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Ab initio modeling of molecular electronics is nowadays routinely performed by combining the 
Density Functional Theory (DFT) and Nonequilibrium Green function (NEGF) techniques. This 
method has its roots in the current formula given by Meir and Wingreen, and we discuss some 
applications and accompanying pitfalls and restrictions of this approach. Quite recently papers have 
begun to appear where inelastic effects are considered, and we illustrate these new developments by 
describing our own work on transport in atomic gold wires. 



I. INTRODUCTION 



A. Some background, and a few historical remarks 

In my previous two reviews in this conference series [1, 2] I have discussed applications of a certain current formula, 
due to Meir and Wingreen[3], to mesoscopic structures. The present review will continue along the same lines. As 
was pointed out in my review in 2002, this formula "is gaining widespread use" ; this has indeed been the case also 
in the intervening three years: more than 350 papers have appeared applying this formula or its variants, and an 
exhaustive review is clearly out of the question. I have therefore restricted myself to two topics which I believe to 
be very important, not only in pure academic sense but also as far as applications are concerned. Specifically, I will 
describe ab initio modeling of molecular electronics, where one combines the density functional theory (DFT) with the 
nonequilibrium Green function technique (NEGF), and, secondly, I address the difficulties, challenges, and successes 
in the very recent efforts where inelastic effects have been included in the modeling. 



B. The basic equations, and their limitations 

The details of the derivation of the Meir- Wingreen formula have been given several times before [4, 5], and there 
is no need to repeat them here, and we just summarize the basic philosophy. In the spirit of the scattering theory 
approach to conductance [6-8], one considers a system (say, a molecule, or a quantum dot) which is coupled via ideal 
leads to reservoirs, which are so large that they can be described by an equilibrium distribution. The nonequilibrium 
aspects are introduced via the following construction, pioneered by Caroli et al.[9-12]: in the infinite past the various 
subsystems are separated, with their respective chemical potentials, and the couplings between the various subsystems 
are turned on adiabatically. The couplings are treated perturbativcly, but to all orders. The double contour of Kcldysh 
enters because only the state in the remote past is known. No assumptions about the smallncss of the couplings are 
needed. One may have to perform a self-consistent calculation of the various parameters describing the system. Below 
we have some comments concerning this construction. 

Let us next consider some generic Hamiltonians: H = Ht, + Hr + Ht + H ccn , or, explicitly: 



H 



^ ek < aC k,a C k,a 
k,a£L/R 



:ncl a d n +h.C. 



k,a£L/ R;n 



H ccn [{dj, {4}] , 



(1) 



where the central part Hamiltonian must be chosen according to the system under consideration. The operators 
{d n }, {d^} refer to a complete set of single-particle states of the central region. The derivation of the basic formula 
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for the current does not require an explicit form for H ccn ; the actual evaluation of the formula of course requires this 
information. We write H ccn = e n dj l d n + H lrit , where H lrA could be electron-phonon interaction, 



or an Anderson impurity: 



tffnV^E^^E^'Waq + S] . (2) 



tf4 = tf£< T d m, T <i<W' ( 3 ) 



or, perhaps, some other model for the Coulomb interaction. 
The current operator for the (say) left lead is 



lL = 'J E [- V kL;ncl L d n + V k * L . n dlc kL \ . (4) 



Many physically relevant observables can be expressed in terms of expectation values of the current operator, or its 
higher powers. For example, one can show[4, 5] that the current leaving the left contact is 

(I L ) = J L (t) = ~J* dt x j glmTrle-^-'^Mx,*) 

x [G<(t,t 1 ) + /«(e)G r (t,t 1 )] } . (5) 

Here the Green functions are defined by 

G< m {tM) = i^faKto) (6) 
G r nm (t,h) = -ieit-h)^),^)]), (7) 

T mn describes the coupling between the central region and the contacts, and /£(e) is the equilibrium distribution 
function of the left contact. The trace in Eq.(5) is over the elements of the product of the T- and G-matrices. 
Further, in Eq.(5) we allow for time-dependent couplings or single-particle energies. Thus, Eq.(5) is applicable for 
charge pumping, which has received a lot of recent attention (sec, e.g., Aono[13] or Kashchcyevs ct al[14]). In the de- 
limit the time-dependence drops out from the r-matrix, and the time-integral generates the energy-representations 
of the G-matrices, which now depend only on the time-difference. After symmetrization (in dc-limit the current 
through left and right barriers are equal; in our notation Jr = —Jl), Eq.(5) then reduces to the result of Meir and 
Wingreen[3]: 

1 = |/|^{[r'(.)-r«( e )]G<( e) 

+ [/EMr'M - f„C)r%)] [G» - &■(<)] } . (8) 

The expressions (5) and (8) are the central formal results whose consequences we explore in this review. They 
are formally exact, and give the tunneling current for an interacting system coupled to noninteracting contacts (or, 
more precisely, for contacts which can be described by an effective single-body Hamiltonian) . The total current 
flowing in the outside circuit may contain contributions from displacement currents, and these must be considered 
separately, see, e.g., Anantram and Datta[15], or Wang et al.[16] . The physical interpretation of (8) can be made 
more transparent by writing it in an alternative form (consider here only current through the left barrier): 

Jl = IJ deTr{S^<( e )G>( £ )-S i '>( e )G<( £ )}, (9) 

where the left self-energy gives the coupling to the left lead. The first term in (9) gives the current flowing through 
the left barrier towards the central region (because it is proportional to G > , i.e., the empty states in the central 
region), while the second term gives the current flowing through the left barrier towards the left contact (because it 
is proportional to G < , i.e., the occupied states in the central region. Likewise, the self-energies are proportional to 
the occupied lead states, and the empty lead states, respectively. 
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It should also be noted that these equations only define the starting point of any calculation: to get physical results 
one must evaluate the correlation function and the retarded/advanced Green function, Eqs.(6) and (7), respectively. 
These functions obey the Keldysh equation, and the (nonequilibrium) Dyson equation: 

G< = (l + G r E r )G^(l + E G°) + G r £ < G a , (10) 
G r = G r + G r Z r G r . (11) 

The first term in the equation for G < represents a boundary term, which is often omitted in discussions of a stationary 
situation, since it only affects the transient behavior. The success of the theory depends on whether one can construct 
a self-energy functional that captures the essential physics, and that a good solution can be found for these coupled 
equations. Both of these steps may be hard indeed. 

While the above approach is widely used, one should not forget its inherent limitations. Consider first the idea 
of partitioning the system into noninteracting leads, and an interacting central region. This separation is crucial for 
obtaining the compact result for the current. But what are the physical criteria for selecting the boundaries between 
the three regions? In our first ab initio calculations of inelastic effects in atomic gold wires [32], to be discussed below, 
we allowed only the four atoms constituting the wire to move around their equilibrium positions and thereby define 
the interacting central region, while in later work [33] several additional atoms from the contact region were also 
taken to be a part of the central region. Including the atoms from the contacts changes the results only slightly, but 
in our case it was the computational effort that dictated how many contact atoms were included, not an a priori 
physical criterion. In a similar vein, allowing the Coulomb interactions to be present only in the central region appears 
somewhat arbitrary. Screening will of course diminish the strength of the Coulomb interactions once the contact has 
broadened into the bulk-like electrode, but once again it is not clear how large the interacting region should be. There 
are also certain subtleties related to overall charge neutrality of the system in nonequilibrium conditions, see Mera 
et. al [17] for some illustrative model calculations. 

The assumption of having noninteracting leads at their respective chemical potentials also hides some important 
physics. The hot charge carriers injected into the collecting electrode are assumed to dissipate their energy, without 
heating up the contacts. There are no terms in the Hamiltonian to describe energy relaxation in the contacts, nor are 
there terms which describe the coupling of the contacts to the surrounding thermal bath. 

The idea of having the system separated into independent parts in infinite past is important in the formal devel- 
opment of the theory: one needs a well-defined initial state which is then modified by external perturbations. In 
an experiment one does not usually introduce the couplings in an adiabatic way at a fixed external bias, and watch 
the build-up of the current flow until it settles to a steady state, rather the entire system is in equilibrium until the 
driving voltage is coupled to the system. This has been recently recognized by Stefanucci and co-authors[18, 19], who 
develop an approach based on time-dependent density functional theory, which avoids the partitioning. These results 
will be discussed elsewhere in these proceedings. 

II. USES AND MISUSES OF THE MEIR-WINGREEN FORMULA 

A. Mean-field theories 

The calculation of the current using (8) together with the noncqilibrium Keldysh and Dyson equations (10) and 
(11), respectively, is a very difficult task for any nontrivial interaction (hidden in the self-energies), and it is natural to 
look for physically motivated simplifications. Let us first consider the case where the interactions are of the mean-field 
type, such as e.g. in the density-functional approach. Mean-field self-energies are one-point functions, with vanishing 
Keldysh components, and the Keldysh equation becomes a formally explicit expression for G < , 

G< = G r i[T L f L + r R f R ]G a . (12) 

Note however that the G r ' a may contain some functional of G < as a parameter (the level occupation is a typical 
example, n n = —iG < (t,t) nn ; another example is Eq.(29) below), and a self-consistent calculation may be needed. 
Nevertheless, Eq.(12) leads to a very compact expression for the current, as we now shall show. Using (12) and the 
relation G r G a = -G r (£ a - S r )G a = iG r TG a (here T = T L +T R ), one finds that the current formulas (8) and 
(9) reduce to 

Jl = IJ deT tot (e)[/°(e)-/°(e)], (13) 

where 



Ttot(e) = Tr{r i (e)G r (e)r fl (e)G a (e)} (14) 
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is the transmission probability. Eq.(13) is used very frequently and occasionally callously: one must bear in mind 
that it only holds for a noninteracting, or a mean-field theory. Nevertheless, it is extremely useful, and forms the 
basis even for some commercial software, which is being used for first-principles modeling of molecular electronics. 



B. Conservation laws 



As emphasized above, the present approach becomes useful the moment one has constructed a self-energy functional. 
But not all self-energy functionals are allowable: they must be such that current conservation is obeyed. This 
important consistency check can be formulated as follows [20]. The sclf-encrgy appearing in the Keldysh equation 
contains contributions both from the coupling to the leads, and from the interactions in the central region. We thus 
write Stot = S; nt + YlaeL/R Using the Keldysh equation repeatedly, one can show that 

Tr{S< t G>-S> t G < } =0. (15) 
Using this condition, one can show that the current conservation condition ^ Q J a = leads to 

J deTr{S< t (e)G>(e)-E? t (e)G<(e)}=0. (16) 

In a practical calculation the sclf-encrgy is known often only numerically and Eq.(16) can be used to check that 
numerical errors have not led to an unphysical result. Using the language of kinetic theory, Eq.(16) says that the 
integrated collision term must vanish - a condition familiar from Boltzmann-like theories. A typical example of a 
"good" self-energy is the self-consistent Born approximation with free equilibrium phonons (which we use below in 
our discussion of inelastic effects in Au-nanowircs): a simple calculation shows that the condition (16) is identically 
satisfied in this case. Calculations considering an interacting and nonequilibrium phonon distribution are scarce, even 
though initial steps towards these goals are being taken[29]. 



C. The wide-band limit 



In the wide-band limit, when the energy-dependence of the T's can be neglected, the Meir-Wingreen formula assumes 
a particularly simple form (for simplicity we consider here a one-state model for the central region): 

J = l^^I J Mh(e) - f R (e)]A(e), (17) 

where A(e) = — 2ImG r (e) is the interacting spectral function. Usually an exact evaluation of the spectral function 
is not possible, and one looks for approximate methods. An often used approach is to start from the atomic limit, 
i.e., an isolated central region, for which the Green function can be calculated, either analytically or perhaps by an 
exact diagonalization, and then broaden the sharp levels by a phcnomcnological width. Important examples include 
a model where a single level is coupled to phonons, for which the Green function is [21] 

G r (t) = -i8(t) cxp[-it(e - A) - $(t)], (18) 

with 



E 



M 2 



M 2 

$ W = V^[A^(l-e i; ^*)+A^(l-e-^*)], (19) 
q ^ 

or an isolated Anderson impurity, for which one has 

G*(e) = {ns) (20) 
e - £cr - U e - e a 

A tempting approach would be to add a factor —Tt/2 in the exponent in the phonon problem, or a factor iY /2 in 
the denominator in the Anderson impurity problem. This would however lead to erroneous physics: in the phonon 
case the heuristic procedure implies that the Fermi seas in the contacts has been ignored (for an extended discussion 
on this issue, see Ref.[30]), and in the Anderson impurity case higher correlations beyond Hartree-Fock level would 
be lost (sec Sect. 12.7 in Rcf.[5]). These simple examples arc mentioned here in order to stress that extreme care is 
needed in treating the interacting problems, and that seemingly innocent steps may have far-reaching consequences. 
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III. DENSITY FUNCTIONAL METHOD FOR NONEQUILIBRIUM ELECTRON TRANSPORT 

In this section I briefly discuss some of the generic features of theories that combine some ab initio electron structure 
theory and a noncquilibrium transport theory In the literature one can find several implementations, here I focus on 
the TRANSIESTA code, developed at the Technical University of Denmark [2 2]. 

Most electronic structure calculations are restricted in the sense that the geometry must be finite, or periodic, and 
that the electronic system is in equilibrium. The situation is very different, say, in the case of molecular electronics: 
here a small subsystem lacking translational invariancc couples to semi-infinite leads and the electronic subsystem 
can be far from equilibrium. Ideally one should describe the whole system (the central region and electrodes) on 
equal footing. As is well-known, the Density Functional Theory gives the exact electronic density and total energy, 
if the exact exchange-correlation functional was known. Since this is not the case, one must resort to approximate 
forms of the functional, such as the local-density approximation (LDA), or the generalized gradient approximation 
(GGA), or something else. There is no theory to say which (approximate) functional is the best, rather the choice is 
made based on painstaking tests, and comparisons in some limits where alternative methods, or experiments, can give 
benchmarks. In an attempt to extend DFT to nonequilibrium situations one must go one step further: the Kohn-Sham 
single-particle wave-functions are used when calculating the current. Thus, genuinely many-particle effects, such as 
the Kondo effect, are excluded from the treatment. Inelastic effects can be included, as discussed in the next section. 
An improvement of the present approach could conceivably be reached by the current-density formalism[23]. 

At the core of the DFT-NEGF implementation described in Ref.[22] is the SIESTA code[24]. This approach has 
been tested in a large number of applications with excellent results, and it has many technical advantages because 
of the employed finite range orbitals for the valence electrons: not only do the numerics get faster but also the 
system partitioning becomes unambiguous. The SIESTA approach can be extended to nonequilibrium by using a 
nonequilibrium electron density as an input; the nonequilibrium electron density is calculated with the help of NEGF: 

n(x) = -iG < (x = x',t = t') = [ ^G < (x = x',e), (21) 

J 2m 

where G < follows directly from the Keldysh equation, because the self-energy is a known function for mean-field 
theories. Explicitly, we recall that the lesser self-energy has the structure 

£< = i(T L f L +T R f R ), (22) 

where the occupation factors of the left and right contacts contain the information about the potential difference. 
Hence, all that one needs are the retarded and advanced Green functions, and these arc obtained by evaluating 

G r ' a (E) = [EI±ir)-H}~ 1 , (23) 

where 



H = V* H c V R , (24) 




with an obvious nomenclature. The semi-infinite left and right leads are accounted for by the self-energies Yi L / R , see, 
e.g., the book by Datta[25]. Importantly, to determine Vjj, Vj,, or He one does not need to evaluate the density 
matrix outside the L — C — R region, if the L — C — R region is defined so large that all screening takes place inside it. 
Summarizing, and somewhat simplifying, the TRANSIESTA iterative loop consists of the steps 

initialn(x) => SIESTA => ^ks(x) => NEGF => newn(x), (25) 

and the iteration is repeated until convergence is achieved for the desired quantity, such as the current for a given 
voltage difference. For a detailed description of many of the technical details suppressed here we refer to the paper 
by Brandbygc et aZ. [22]. The scheme outlined above has been applied by a large group of researchers to many 
specific physical systems. Occasionally the agreement with experiments reaches a quantitative level, which is indeed 
very satisfying, while sometimes the predicted current can be orders of magnitude too large. At the moment of 
writing this review, there is no consensus of whether the discrepancies are due to poorly controlled experiments, bad 
implementations of the DFT-NEGF scheme, or due to an inadequacy of the entire concept. A possible cause for 
the discrepancy has very recently been identified by Toher et al.[34], who suggest that self-interaction corrections 
(which are not included in the GGA-LDA underlying most theoretical work, including ours) could remedy some of 
the problems. Nevertheless, a lot of research remains to be done. 
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IV. INELASTIC SCATTERING AND LOCAL HEATING IN ATOMIC GOLD WIRES 

The issue of vibrational effects in molecular electronics has recently drawn a lot of interest because inelastic scat- 
tering and energy dissipation inside atomic-scale conductors are of paramount importance for device characteristics, 
working conditions, and their stability [26-28, 30]. Inelastic effects arc important, not only because of their poten- 
tially detrimental influence on device functioning, but also because they can open up new possibilities and operating 
modes. Vibrational effects are often visible in the measured conductances of nanoscopic objects; here we focus on 
recent experimental studies on free standing atomic gold wires [31]. Agrait and co-workers used a cryogenic STM tip 
to first create an atomic-scale gold wire (lengths up to seven gold atoms have been achieved), and then measured 
its conductance as a function of the displacement of tip, and the applied voltage. The data showed clear drops of 
conductance at a certain voltage, and the interpretation was that an excitation of an inelastic mode was taking place, 
leading to enhanced back-scattering, and hence drop in the conductance. It should be pointed out that opening a 
new vibrational mode in the atomic scale conductor does not necessarily lead to a decrease in conductance (one can 
envisage various assisted processes) , and a proper theory should be able to predict conductance enhancement as well, 
whenever the physics dictates so. 

Here I will briefly describe our recent work[32] on inelastic effects in the kind of wires studied by Agrait et al.[31]. 
An important aspect of our work is that we go beyond lowest order perturbation theory in the electron-vibration 
coupling, and therefore polaronic effects can be included. We also address the issue of phonon heating, albeit within 
a phenomenological model, as discussed below. 

The calculational method consists of three steps, (i) The mechanical normal modes and frequencies of the gold 
chain are evaluated, (ii) The electronic structure and electron- vibration coupling elements arc evaluated in a localized 
atomic-orbital basis set. (iii) The inelastic transport is evaluated using NEGF by using a self-consistent Born approx- 
imation self-energy in the Dyson and Keldysh equations for the respective Green functions. The electrical current 
and the power transfer are then evaluated with (here, for the left lead; for a detailed derivation, see Ref.[20]) 

h = ~f det L {e) 

t L (c) = lt{Sj( e )G>( e )-Ej(e)G<( £ )}, 
where Hartrec and Fock parts of self-energy components are 

S i?,r = { y 2 I ^ M A Tr[G<(e / )M Aj 
A J 

E ff '< = 
sF ' r W = *E/£ MA [*-OG<(e') 

+D r Q (e e')G r (e') + D<(e - e')G r (e')]M 
SF ^) = i£/^M^ <(u;-u/)G<(e')M\ 

A J 

Here the vibrational modes are labeled by A, and fl\ is the corresponding eigcnfrcqucncy. It is worth noting that 
the lack of translational invariance makes the retarded Hartree term non-zero, and potentially important. Also, at 
this stage the phonon propagators are undamped - an approximation that merits further investigation. The coupled 
equations are iterated until convergence is achieved, and in the following we give some representative results. Let us 
consider a linear four-atom gold wire under two different states of strain, as shown in Fig. 1. We calculate the phonon 
signal in the nonlinear differential conductance vs. bias voltage for two extremal cases: the energy transferred from 
the electrons to the vibrations is either (i) instantaneously absorbed into an external heat bath, or (ii) accumulated 
and only allowed to leak via electron-hole pair excitations. These limits are referred to as the externally damped 
and externally undamped cases, respectively. The next two figures show the corresponding results. Since a typical 
experiment is done at low temperatures, the mode occupation in the externally damped case vanishes, N\ ss 0. In the 
externally undamped case the mode occupation N\ is an unknown parameter entering the electron-phonon self-energy, 
and additional physical input is necessary to determine this parameter. We argue as follows. Since the system is in a 
steady state, the net power transferred from the electrons to the device must vanish, i.e., Pj, + Pr = 0. Using Eq.(27) 
one then obtains the required constraint on N\. This procedure works in a straightforward way if there is only a single 



(26) 

(27) 
(28) 



(29) 
(30) 

(31) 
(32) 
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FIG. 1: Geometry of a four-atom gold wire under two different states of stress. The dominating alternating bond-length modes, 
which cause the inelastic scattering, are shown schematically with arrows. For the shorter wire only one mode is active, while 
the elongated wire has two active r~ ~ J ~~ /T1 J J £ ' ™ ~ r ror,lx 




FIG. 2: Differential conductance and its derivative for the four-atom gold wire at two different tensions in the case where the 
oscillators are externally damped. Reproduced from Ref. [32]. 



active mode, but if several modes are present, a more detailed theory of how the phonon modes equilibriate would 
be needed. The bottom panel of Fig. 3 shows how the occupation of the leading phonon mode, or, equivalently, the 
effective temperature, changes as a function of bias. When comparing to the experiments of Agrait et al.[31], one sees 
that the externally undamped model is in near quantitative agreement with the data: the conductance drop at the 
onset of inelastic scattering, and the slope after the drop are very well reproduced. We view this as strong evidence of 
the presence of heating in the experiment, but at the same time recognize the need for a detailed microscopic theory 
including phonon-phonon interactions. 



V. CONCLUSIONS 



I have reviewed some of the post-2002 developments in applying NEGF to modelling of transport in mesoscopic 
systems. The common theme in my review has been to focus on "real devices" , which may have "real applications" . 
I find it very pleasing that the NEGF technique, often regarded as an academic exercise most suited for theoretical 
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FIG. 3: (a) Differential conductance and its derivative for the four-atom gold wire at two different tensions in the externally 
undamped limit. Only the most important mode is included in this calculation, (b) Mode occupation vs. bias voltage. 
Reproduced from Ref.[32]. 



games, is now becoming a strong tool in the analysis of practical devices, even in industrial context. At the same 
time there is still much room for theoretical refinements, and I'm convinced that in the coming years we will witness 
significant progress in this field, both abstract and practical. 
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